Comparison distance-based versus model-based approach

Packages

Custom functions

#--------------------------------------------------
# Plot ordinations from rda(), cca(), or cmdscale()
#--------------------------------------------------

plot_ord <- function(ord,comm, reverse_axis2 = FALSE){
  
  site_scores <- scores(ord, display = "sites") %>%
    as.data.frame() %>%
    set_names("Axis1", "Axis2")
  
  if(reverse_axis2){
    site_scores$Axis2 <- -site_scores$Axis2 
  }
  
  if(any(class(ord) %in% c("rda","cca"))){
    
   sp_scores <- scores(ord, display = "sp") %>%
     as.data.frame() %>%
    set_names("Axis1", "Axis2")
    
  }else{# PCoA
    
    sp_scores <- wascores(site_scores, comm) %>%
      as.data.frame()
  }
 
  ggplot() +
    geom_point(data = site_scores, aes(x = Axis1, y = Axis2)) +
    geom_path(data = site_scores, aes(x = Axis1, y = Axis2)) +
    geom_text_repel(data = site_scores, aes(x = Axis1, y = Axis2, label = rownames(site_scores)), nudge_y = -0.1) +
    geom_segment(data = sp_scores, aes(x = 0, y = 0, xend = Axis1, yend = Axis2), col = "grey") +
    geom_text_repel(data = sp_scores, aes(x = Axis1, y = Axis2, label = rownames(sp_scores)), col = "grey") +
    coord_fixed(ratio = 1) +
    theme_minimal()
   
}


#--------------------------------------------------
# Plot ordinations from boral
#--------------------------------------------------

plot_boral <- function(boral_mod){
  
  # Retrieve data to do the same plots as before
  ord_data <- gg_lvsplot(boral_mod) %>%
    layer_data(., 2) %>%
    select(x, y, label) %>%
    mutate(group = if_else(str_detect(label, "sp"), "species", "sites"))
  
   ggplot() +
    geom_point(data = ord_data %>% filter(group == "sites"), aes(x = x, y = y)) +
    geom_path(data = ord_data %>% filter(group == "sites"), aes(x = x, y = y)) +
    geom_text_repel(data = ord_data %>% filter(group == "sites"), aes(x = x, y = y, label = label), nudge_y = -0.1) +
    geom_segment(data = ord_data %>% filter(group == "species"), aes(x = 0, y = 0, xend = x, yend = y), col = "grey") +
    geom_text_repel(data = ord_data %>% filter(group == "species"), aes(x = x, y = y, label = label), col = "grey") +
    coord_fixed(ratio = 1) +
    theme_minimal() +
     xlab("Latent variable 1") +
     ylab("Latent variable 2")
   
}

Legendre et Gallagher 2001

Reproduce Fig. 4

Using decostand :

  • profile transformation: Y.tr = decostand(Y, “total”)
  • chord transformation: Y.tr = decostand(Y, “norm”)
  • log-chord transformation: Y.tr = decostand(log1p(Y), “norm”)
  • Hellinger transformation: Y.tr = decostand(Y, “hellinger”)
  • chi-square transformation: Y.tr = decostand(Y, “chi.sq”)

dist.ldc can also be used

\s

Figure 2:

  • Add variance explained to axes labels
  • Do log-chord

Comparison of Fig. 4 with JSDM

Comparison with boral

Normal

# Normal - without row effects
#----------------------------

boral_normal <- boral(y = comm, family = "normal", lv.control = list(num.lv = 2), row.eff = "none")

boral_p1 <- plot_boral(boral_normal)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Normal - without row effects")

# Normal - fixed row effects
#----------------------------

boral_normal_fixed <- boral(y = comm, family = "normal", lv.control = list(num.lv = 2), row.eff = "fixed")

boral_p2 <- plot_boral(boral_normal_fixed)+
  scale_y_reverse() +
  ggtitle("Normal - fixed row effects")

# Normal - random row effects
#----------------------------

boral_normal_rand <- boral(y = comm, family = "normal", lv.control = list(num.lv = 2), row.eff = "random")

boral_p3 <- plot_boral(boral_normal_rand)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Normal - random row effects")

# Hellinger - Normal - without row effects
#----------------------------
comm_tr <- decostand(comm,"hellinger")

boral_normal_hell <- boral(y = comm_tr, family = "normal", lv.control = list(num.lv = 2), row.eff = "none")

boral_p4 <- plot_boral(boral_normal_hell)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Hellinger-transformed - Normal - without row effects")

# Normal - fixed row effects
#----------------------------

boral_normal_hell_fixed <- boral(y = comm_tr, family = "normal", lv.control = list(num.lv = 2), row.eff = "fixed")

boral_p5 <- plot_boral(boral_normal_hell_fixed)+
  scale_y_reverse() +
  ggtitle("Hellinger-transformed - Normal - fixed row effects")

# normal - random row effects
#----------------------------

boral_normal_hell_rand <- boral(y = comm_tr, family = "normal", lv.control = list(num.lv = 2), row.eff = "random")

boral_p6 <- plot_boral(boral_normal_hell_rand)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Hellinger-transformed - Normal - random row effects")
\s

Figure 3:

\s

Figure 4:

Poisson & Negative binomial

# Poisson - without row effects
#----------------------------

boral_poisson <- boral(y = comm, family = "poisson", lv.control = list(num.lv = 2), row.eff = "none")

boral_p1 <- plot_boral(boral_poisson)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Poisson - without row effects")

# Poisson - fixed row effects
#----------------------------

boral_poisson_fixed <- boral(y = comm, family = "poisson", lv.control = list(num.lv = 2), row.eff = "fixed")

boral_p2 <- plot_boral(boral_poisson_fixed)+
  scale_y_reverse() +
  ggtitle("Poisson - fixed row effects")

# Poisson - random row effects
#----------------------------

boral_poisson_rand <- boral(y = comm, family = "poisson", lv.control = list(num.lv = 2), row.eff = "random")

boral_p3 <- plot_boral(boral_poisson_rand)+
  scale_y_reverse() +
  scale_x_reverse() +
  ggtitle("Poisson - random row effects")

# Negative binomial - without row effects
#--------------------------------------

boral_negbin <- boral(y = comm, family = "negative.binomial", lv.control = list(num.lv = 2), row.eff = "none")

boral_p4 <- plot_boral(boral_negbin) +
  scale_y_reverse() +
  ggtitle("Negative binomial - without row effects")

# Negative binomial - fixed row effects
#--------------------------------------

boral_negbin_fixed <- boral(y = comm, family = "negative.binomial", lv.control = list(num.lv = 2), row.eff = "fixed")

boral_p5 <- plot_boral(boral_negbin_fixed) +
  scale_y_reverse() +
  ggtitle("Negative binomial - fixed row effects")

# Negative binomial - random row effects
#--------------------------------------

boral_negbin_rand<- boral(y = comm, family = "negative.binomial", lv.control = list(num.lv = 2), row.eff = "random")

boral_p6 <- plot_boral(boral_negbin_rand) +
  scale_y_reverse() +
  ggtitle("Negative binomial - random row effects")
\s

Figure 5:

Residuals of the models

  • Normal - without row effects
## NULL
\s

Figure 6:

\s

Figure 7:

\s

Figure 8:

\s

Figure 9:

  • Poisson - without row effects
## NULL
\s

Figure 10:

\s

Figure 11:

\s

Figure 12:

\s

Figure 13:

  • Negative binomial - without row effects
## NULL
\s

Figure 14:

\s

Figure 15:

\s

Figure 16:

\s

Figure 17:

Comparison with mvabund

\s

Figure 18:

Session info

## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.15.5
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] boral_1.8       coda_0.19-3     mvabund_4.0.1   vegan_2.5-6    
##  [5] lattice_0.20-35 permute_0.9-5   ggboral_0.1.7   patchwork_0.0.1
##  [9] ggrepel_0.8.2   forcats_0.3.0   stringr_1.4.0   dplyr_1.0.2    
## [13] purrr_0.3.2     readr_1.1.1     tidyr_1.0.0     tibble_2.1.3   
## [17] ggplot2_3.3.0   tidyverse_1.2.1
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-137     lubridate_1.7.4  httr_1.4.1       tools_3.5.0     
##  [5] backports_1.1.4  R6_2.4.1         R2WinBUGS_2.1-21 mgcv_1.8-23     
##  [9] questionr_0.7.0  colorspace_1.4-1 withr_2.1.2      tidyselect_1.1.0
## [13] compiler_3.5.0   cli_2.0.2        rvest_0.3.4      xml2_1.2.2      
## [17] labeling_0.3     bookdown_0.13    scales_1.0.0     mvtnorm_1.1-0   
## [21] digest_0.6.25    rmarkdown_1.15   pkgconfig_2.0.3  htmltools_0.3.6 
## [25] highr_0.7        rlang_0.4.8      readxl_1.1.0     rstudioapi_0.9.0
## [29] shiny_1.3.2      generics_0.0.2   jsonlite_1.6     magrittr_1.5    
## [33] Matrix_1.2-14    Rcpp_1.0.5       munsell_0.5.0    fansi_0.4.1     
## [37] abind_1.4-5      lifecycle_0.2.0  stringi_1.4.6    yaml_2.2.0      
## [41] MASS_7.3-51.4    plyr_1.8.4       fishMod_0.29     grid_3.5.0      
## [45] parallel_3.5.0   promises_1.0.1   crayon_1.3.4     miniUI_0.1.1.1  
## [49] haven_2.1.1      hms_0.4.2        knitr_1.24       pillar_1.4.2    
## [53] boot_1.3-20      corpcor_1.6.9    reshape2_1.4.3   codetools_0.2-15
## [57] R2jags_0.5-7     glue_1.4.2       evaluate_0.14    modelr_0.1.5    
## [61] vctrs_0.3.5      rmdformats_0.3.5 httpuv_1.5.2     cellranger_1.1.0
## [65] gtable_0.3.0     assertthat_0.2.1 xfun_0.9         mime_0.7        
## [69] tweedie_2.3.2    xtable_1.8-4     broom_0.5.2      later_0.8.0     
## [73] rjags_4-10       cluster_2.0.7-1  statmod_1.4.34   ellipsis_0.3.0

Aurélien Boyé

18 décembre, 2020